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Networks — collections of interacting elements or nodes — abound in the natural and manmade worlds. For 
many networks, complex spatiotemporal dynamics stem from patterns of physical interactions unknown to us. 
To infer these interactions, it is common to include edges between those nodes whose time series exhibit suffi- 
cient functional connectivity, typically defined as a measure of coupling exceeding a pre-determined threshold. 
However, when uncertainty exists in the original network measurements, uncertainty in the inferred network 
is likely, and hence a statistical propagation-of-error is needed. In this manuscript, we describe a principled 
and systematic procedure for the inference of functional connectivity networks from multivariate time series 
data. Our procedure yields as output both the inferred network and a quantification of uncertainty of the most 
fundamental interest: uncertainty in the number of edges. To illustrate this approach, we apply our procedure 
to simulated data and electrocorticogram data recorded from a human subject during an epileptic seizure. We 
demonstrate that the procedure is accurate and robust in both the determination of edges and the reporting of 
uncertainty associated with that determination. 



I. INTRODUCTION 

Many examples of natural and fabricated networks exist in 
the world, including airline networks, computer networks, and 
neural networks. To define a network is, in principle, straight- 
forward: we simply identify a collection of nodes and edges 
CI QUI]- A node (or vertex) represents a participant or actor in 
a network, while an edge represents a link or association be- 
tween two nodes (Fig.[TJ). For example, in an airline network, 
individual airports constitute nodes and direct connections be- 
tween airports identify edges. In a neural network, individual 
neurons and the physical connections between neurons deter- 
mine the nodes and edges of the network, respectively. Hav- 
ing established network representations of these complex sys- 
tems, we may then address pertinent issues, such as the world- 
wide spread of infectious disease through the airline network 
J3l or the effect of cortical lesions on brain dynamics JH] . 

The decision to link two nodes with an edge varies in dif- 
ficulty. In some cases, a known physical connection exists 
between two nodes, and the choice to include an edge is then 
obvious. Does an airline connect two cities or not |@]? Do 
two actors collaborate on a film or not 0,111]? Does a phy sical 
connection exist between two brain regions or not lloL 1 1 Oh ? In 
these cases, the decision to include a link between nodes is 
simple and based on the known association or physical con- 
nection between two nodes. 

In other cases, the interactions between nodes are obscure. 
For example, we may only observe the dynamic activity at in- 
dividual nodes and have no access to the physical connections 
between nodes. In these cases, we may apply coupling mea- 
sures to multivariate time series data associated with the node 
dynamics and attempt to infer their associations, without ex- 
plicit knowledge of their structural connections fTH [l2l [l3ll . 
This approach has proven useful in, for example, climate stud- 
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ies lfl4l [l5ll and also human brain studies, in which the struc- 
tural connections between brain regions remain difficult to 
classify (although perhaps not for long lfl6ll ). 

Viewed from a statistical perspective, two key challenges 
are inherent in this task of network inference: (i) appropri- 
ate interpretation of the coupling results in declaring network 
edges, and (ii) accurate quantification of the uncertainty asso- 
ciated with the resulting network. The simplest - and, indeed, 
most common - method to interpret the coupling results and 
declare network edges involves comparison of the coupling 
strength to a threshold value H El 111 El IS 111 111 111. If 
the coupling strength between two nodes exceeds this thresh- 
old, then we connect these nodes with an edge; otherwise, we 
leave the nodes unconnected. The number of edges in the re- 
sulting network depends critically on the choice of coupling 
threshold (as we illustrate schematically in Fig. Q]). Further- 
more, for a given choice of threshold, we expect a certain 
rate of error in (mis)declaring the connectivity status between 
pairs of nodes. This network uncertainty — intimately tied to 
choice of threshold — is often overlooked. 

How do we choose such a threshold? One strategy is to ap- 
ply a variety of different thresholds and examine the resulting 
collection of networks for robustness as a function of thresh- 
old HH I2H E3I1 . This procedure of redundant anal ysis — 
which some of these authors have recently employed M22I1 — 
is both time consuming and unsatisfying. Instead, a thresh- 
old should be chosen in a principled way, for example one 
that links the choice of threshold with the achievement of a 
pre-specified level of network uncertainty. Such is the goal of 
this paper. In referring to "network uncertainty" many aspects 
of the network structure might be of interest (e.g., connectiv- 
ity, degree distribution, or clustering). Here we will focus on 
the most basic aspect of network confidence: the presence or 
absence of network edges. Our particular goal is to equip the 
process of threshold-based inference of a network with a num- 
ber quantifying the expected rate of falsely declared network 
edges. This number serves as a natural measure of network 
confidence. 

In what follows, we adopt a statistical hypothesis testing 
paradigm to analyze multivariate time series data and create 
network representations of functional connectivity. The gen- 
eral paradigm involves three steps: 1) calculate the strength of 
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coupling between time series data recorded at node pairs, 2) 
threshold each coupling measure through the use of a formal 
statistical hypothesis test, and 3) control the rate of falsely 
declared edges through the use of statistical multiple testing 
procedures. In Section [II] we present a high-level outline of 
this general protocol, while in Section[III] we develop the pro- 
cedure in detail, making specific choices of methodology for 
each step. We apply the protocol to three data sets in Section 
HVI and show, in particular, that appropriate handling of the 
significance tests is vital. Fundamentally, the proposed proto- 
col is a way of constructing functional networks that, rather 
than emerging as the result of some arbitrarily chosen cou- 
pling threshold, are composed of edges selected to achieve a 
guaranteed level of overall network accuracy. That is, it is a 
way of constructing networks with confidence. 




Threshold = 0.05 Threshold = 0.65 Threshold = 0.95 



FIG. 1: The number of edges in this 4-node network depends upon 
our choice of coupling threshold. Each node (gray circle) represents 
a spatial location from which we record time series data. Applying 
a coupling measure to the dynamic activity recorded at node pairs, 
we obtain the normalized (between and 1) coupling values shown 
in (a). If we choose the coupling threshold too low (0.05) we in- 
clude edges between all nodes as in (b). If we choose the coupling 
threshold too high (0.95) we obtain no edges in the network as in 
(d). An intermediate choice of coupling threshold (0.65 in c) yields 
a different network. 



II. GENERAL PARADIGM 

In this paper we are interested in the inference of net- 
works (or, more precisely, graphs) G = (V,E) in which edges 
{i,j} 6 E indicate a coupling (perhaps at nonzero lag) be- 
tween time series Xj[t] and Xj[t] observed at N nodes i,j G V. 
Our primary motivation is the desire to infer networks reflect- 
ing the functional (as opposed to structural) topology of neural 
systems. This goal is reflected in our terminology, as well as 
in the numerical illustrations we present in Section [IV] How- 
ever, the methods we propose — and the underlying principles 
upon which we base the methods — have quite general appli- 
cability. 

Network inference problems come in many varieties. See 
Chapter 7 of l24ll for a recent overview of this highly active 
area. The type of networks we wish to infer are commonly 
called association networks. Broadly speaking, most methods 
proposed to date for the inference of such networks assume 
independent measurements at each node. A primary example 
of this paradigm is the popular problem of inferring Gaus- 
sian graphical models. Methods for doing so include clas- 
sical methods of maximum-likelihood-based testing (e.g., see 
Chapter 6 of l25ll or Chapter 5 of ll26ll ). and more recent meth- 



ods based on multiple tes ting (e .g., l27l 1281 l29ll ) and sparse 
statistical inference (e.g., 1301, 13111 ). 

However, work of this sort invariably assumes independent 
measurements in time rather than temporally correlated time 
series of interest here. And furthermore, most of these meth- 
ods (e.g., classical and those based on sparse inference prin- 
ciples) are not aimed at providing a quantification of uncer- 
tainty with the inferred network. Alternatively, there is also 
a sub-literature on the inference of association networks from 
temporally correlated data (e.g., S [H [H [H). But the 
quantification of network uncertainty does not seem to have 
received much attention there. 

We implement a procedure to create functional topologies 
from multivariate time series data that involves three general 
steps. First, a coupling measure — here, the cross correlation 
— is specified and applied to the data, yielding a noisy indi- 
cation of the functional connectivity between all nodes pairs. 
In the neurological data described below, this measure cap- 
tures the extent of interactions between activity recorded si- 
multaneously at separate spatial locations of the brain. Sec- 
ond, we develop significance tests appropriate for our choice 
of coupling measure, and associate a statistical /?-value with 
each coupling result. Third, we analyze the resulting p-values 
using principles of statistical multiple testing to construct a 
network representation of the functional connectivity. In the 
course of this last step, we determine a number controlling the 
proportion of falsely inferred edges. We present this number 
as a basic and natural measure of network uncertainty. 

Of course, the first step above implements a version of 
the standard approach to constructing such functional net- 
works. In neuroscience, for example, investigators (including 
some of these authors) typically specify a measure of coupling 
and then assign edges between node pairs whose coupling is 
judged to be sufficiently strong. However, determination of 
just how strong is strong enough is invariably ad hoc or, at 
best, driven by "expert judgement". As a result, there is no 
way to annotate the resulting networks with any indication 
of their inherent (in)accuracy. The subsequent steps in the 
proposed approach, therefore, are critically important to pro- 
duce networks accompanied by accurate characterizations of 
their uncertainty. Put another way, we are interested here in 
the propagation of uncertainty in network inference, from the 
original time series data x, [f] to a final assessment of network 
uncertainty. Our statistical hypothesis testing procedure, de- 
scribed above in three steps, achieves this goal. Furthermore, 
our numerical results indicate that it in fact does so in a robust 
fashion. 

We achieve our goal primarily through careful attention to 
the interdependency among each of the three steps. In so do- 
ing, we also demonstrate how lack of such attention can lead 
to nonsensical network uncertainty statements. For example, 
the choice of coupling measure in the first step affects the hy- 
potheses tested in the second step (i.e., the null hypothesis Hq : 
No Coupling, versus Hi : Coupling). The declaration of either 
edge or non-edge for each pair of nodes i, j G V corresponds 
to either rejection of the null hypothesis or a failure to do so, 
respectively. If rejection is determined by comparison of the 
observed coupling values to a threshold, clearly the choice of 
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threshold will affect the network results. But if we also wish 
to propagate uncertainty — from the original time series data, 
through the testing procedure, to the final network inferred - 
it is necessary to construct accurate probabilistic statements 
appropriate for the particular coupling measure we choose. 

In constructing functional networks, we must consider the 
collection of individual hypothesis tests as a whole. We note 
that the classical approach to calibrating individual hypothe- 
sis tests is not appropriate here. If our hypothesis tests are 
conducted at some significance level a then, for each pair 
of nodes i,j G V, we expect an edge to be declared falsely 
between them with probability a. However, since there are 
N(N — l)/2 such tests to be conducted (assuming an undi- 
rected network of nodes), we actually expect a[N(N— l)/2] 
edges to be declared falsely over the network as a whole, as- 
suming independence of tests. This suggests that as the net- 
work size increases we must decrease a to limit the total num- 
ber of falsely declared edges. But this strategy in turn has 
the undesirable effect of decreasing the statistical power with 
which we can detect edges. This is the so-called "multiple 
testing problem" in statistics. 

Alternatively, therefore, we instead focus upon controlling 
the rate of falsely declared edges. Conditional on at least one 
edge being declared, the expected proportion of falsely de- 
clared edges here is equivalent to what is called the false dis- 
covery rate (FDR) in the statistical literature. The control of 
FDR in multiple testing situations, ranging from signal and 
image processing to genome-wide association testing, has be- 
come a de facto standard technique for addressing the multi- 
ple testing problem (see 13211 . for example). In fact, the use 
of FDR control occurs with increasing frequency in the net- 
work literature as well (e.g., l33Tl34ll ). However, there is little 
evidence in this literature that the rates quoted are necessar- 
ily being achieved. We show later that it is quite easy, using 
a seemingly reasonable significance test, to end up with rates 
that are completely unrepresentative. 

III. IMPLEMENTATION OF THE GENERAL PARADIGM 

We described in the previous section three general steps 
to create a functional topology from multivariate time series 
data. In this section, we 1) define our coupling measure, the 
cross correlation, 2) develop appropriate significance tests, 
and 3) integrate these with a technique to account for mul- 
tiple significance tests. In Section|IV]we apply these specific 
protocols to simulated and observed data, and show that the 
choice of the significance test is critical. 

A. Step 1: choice of coupling measure 

The choice of coupling measure between pairs of time se- 
ries permits many alternatives ll3~5ll . We may select a sim- 
ple measure of linear coupling (e.g., the cross correlation 
111 d |H [H H3l or the coherence HI [M HI) or more 
sophisticated coupling measures (e.g., synchronization likeli- 
hood Qjj] , wavelet coherence |4lll . or Granger causality and 



the related directed transfer function J42tl ). In this manuscript, 
we choose to focus on a simple measure of coupling based on 
the cross correlation. Although the general statistical hypoth- 
esis testing paradigm we adopt here can in principle be ap- 
plied to any choice of coupling measure, more sophisticated 
coupling measures may not easily allow for the derivation of 
computationally tractable significance testing procedures. 

Specifically, for a pair of time series Xi[t] and xj[t], the cross 
correlation at lag T is defined as 

CQj [T] = - 1 X> [t] " *■) (Xj [t + X] - Xj) , ( 1 ) 

OiOj(n -2x) fr[ 

where I, and xj are the averages, and a,- and Oj are the stan- 
dard deviations of xj[t] and xj[t), respectively, and n is the time 
series length. This quantity can be computed efficiently over 
a range of lags using Fourier transform methods for convolu- 
tions. In our applications, we first transform the time series at 
each network node to have zero mean and unit variance, af- 
ter which we compute the Fourier transforms x,[co] and Xj[(o], 
multiply the first by the complex conjugate of the second, and 
take the inverse Fourier transform of the resulting product. For 
all of the data considered below, we compute the cross corre- 
lation for x ranging over indices between —100 and 100 (the 
range of x in milliseconds depends on the sampling rate of the 
data, as we describe below). 

Our formal measure of coupling will be the maximal cross 
correlation i.e., Sjj = max x |CCy[x]|, the maximum of the ab- 
solute value of CCjj [x] over x. This measure will serve as our 
statistic for testing whether or not to assign an edge between 
nodes i and j, for each such pair of nodes. 

B. Step 2: Significance test 

Having chosen the test statistic sij, the maximal cross cor- 
relation between jc,-[f] and xj[t], do we include a network edge 
between nodes ; and jl To answer this, we will use su to test 
the null hypothesis that x, [t] and xj [t] are uncorrected (i.e., no 
coupling) against the alternative that they are correlated (i.e., 
coupling). Rather than focusing on testing at a pre-assigned 
significance level, we will instead concentrate first on comput- 
ing an appropriate p-value for each edge. Accurate evaluation 
of the p-values is critical to successful use of the false dis- 
covery rate principles we employ for the network inference 
problem here, as we discuss in Section ITlI CI We compute the 
p-value in two different ways that make different assumptions 
about the coupling results. The first method is an analytic 
measure and specifically designed for our choice of coupling 
measure. The second is more general but computationally ex- 
pensive. We define the measures below and, in the next sec- 
tion, apply each to simulated and observed time series data. 

1. Analytic method 

In this section we propose an analytic method. Frequently 
such methods involve comparison of a test statistic to a nor- 
mal distribution. Following this approach, we would scale sn 
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by an estimated variance, and then compare this scaled quan- 
tity to the standard normal distribution (i.e., with mean and 
variance 1) to calculate a p-value. Here, however, this would 
be naive. 

More specifically, under the stated null hypothesis of no 
coupling the statistic sy should have mean zero. A reasonable 
estimate of the variance of CCy [I] under the null, motivated 
by a result of Bartlett JHHl], is given by 

var(/) = — £ CC H [T]CCjj[x], (2) 

" T=-n 

where the CC^[x] are the autocorrelations of time series k at 
lag X. This estimate takes a non-trivial form because the cross 
correlation will depend on the statistical properties of the un- 
derlying time series, and in particular on the autocovariance 
structure. Spurious cross correlations between the two times 
series are expected even if they are uncoupled l45ll . and this 
variance formula accounts for that. 

Intuitively we might think to use (f2]) to define zy = 

/ \Jvai(ljj), where Zy is the lag corresponding to sij (i.e., 
the lag at which the maximum of the absolute value of the 
cross correlation occurs) and test the significance of the value 
Zij by comparing it to the standard normal distribution. Un- 
fortunately, although standardizing CCy[/] by the estimated 
variance in (O is sensible for any fixed I, use of the standard 
normal distribution with Zij is not appropriate here, as we ex- 
plain and illustrate below. 

Two potential problems exist in using this naive method to 
determine the significance of sy. First, the distribution of the 
cross correlation Cy[x] - strictly speaking - is normal only 
in the asymptotic case of large sample size n. In finite sam- 
ples this approximation can be poor, particularly since the 
cross correlations are bounded between -1 and 1 while the nor- 
mal distribution varies over an unbounded range. Second, we 
choose sn as an extremum of the cross correlation; therefore, 
we must account for this choice when testing the significance 
of this statistic. That is, even in cases where the distributions 
of the cross correlations Cy [x] are well-approximated by the 
normal distribution, their extrema will not be normally dis- 
tributed, and so p-values calculated using this distribution will 
be inaccurate. 

To address both of these issues, we propose a more appro- 
priate analytic method: the extremum method. We start by 
applying the well-known Fisher transformation H46I1 to each 
CQj[x], yielding 

which should more closely follow a normal distribution than 
the original CCy[x]. Since this transformation is monotone 
and symmetric about zero, the lag Zy maximizing |CCy[x] | will 
also maximize |FCCy [x] |. Let sfj be the Fisher transformation 
of s^, which we propose to use instead of sy. 

Next, we use results from extreme value theory to approx- 
imate the distribution of our new test statistic. We scale the 
values FCCij [x] over x by their empirical standard deviations 



{var{FCCij)) 1 / 2 so that the resulting scaled values should ap- 
proximately follow a standard normal distribution. If there 
were no dependency within the time series x, [f], and instead 
we observed i.i.d. sequences at each node ;, then the appro- 
priate standard deviation is known to be (n — 3)~'/ 2 l46ll . But 
given the dependency in our time series data, we expect that 
the true standard deviations may differ from this value, and so 
we choose to estimate them empirically. 

The scaled value zfj = 4j/(var(FCCij)) 1 / 2 can be expected 
to behave like the maximum of the absolute values of a se- 
quence of standard normal random variables. Using estab- 
lished results for statistics of this form, we obtain therefore 
that 

P[z]«exp(-2exp(-a„(z -£„))) , (4) 

where P[z] = Pr{zf,- < z}, a„ = y/2 log n and b„ = a n — 
(2a n )~ l (loglogn +log47t). A derivation of (0), which holds 
in the asymptotic sense of large n, is provided in Appendix 1 . 
For the case of « = 201, as in all of our numerical results 
below, a„ = 3.2568 and b n = 2.6121. Using the approxima- 
tion above, it is straightforward to calculate p-values for the 
rescaled test statistics zfj- 

Intuitively, the extremum method accounts for our choice of 
a maximum cross correlation. By virtue of the Fisher transfor- 
mation, the values FCCy[x] will be approximately normally 
distributed. But because we have chosen sfj as the maximum 
of the absolute value of the FCCs, we expect its value to be 
skewed towards the tail of the normal distribution. If we had 
chosen instead any other lag than that maximizing the cross 
correlation, then the corresponding value CCy (and hence 
FCCij) would be smaller. Therefore, our definition of sfj pro- 
duces /^-values that, if computed from the normal distribution, 
are biased in the sense of being inappropriately small. From 
the perspective of our network inference task, this means that 

— for any given choice of threshold — we will be more lib- 
eral in our assignment of edges than we should be. The dis- 
tribution in (0| essentially corrects for this bias, by explicitly 
accounting for our use of the maximum. 

2. Frequency domain bootstrap method 

The previous method provides an analytic formula for test- 
ing the significance of sy. In utilizing this formula, we make 
specific asymptotic distributional assumptions about the test 
statistic — the maximal cross correlation. These assumptions 
are likely to be only approximate idealizations of the correla- 
tion results emerging from, for example, a complicated phys- 
ical system like the human brain. A method to test the signif- 
icance of sn that requires fewer assumptions is desirable. The 
final method we introduce — the frequency domain bootstrap 

— satisfies this desire, but is computationally expensive. 

As the name indicates, the method consists of applying the 
bootstrap principle (e.g., PiTll ). but in the spectral domain; 
methods of this sort were first proposed in 14811 . We calcu- 
late our frequency domain bootstrap through the following 
steps. First, we compute the power spectrum (Hanning ta- 
pered) of each time series in the network. We then average 
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these power spectra from all time series, and smooth the re- 
sulting average spectrum (moving average of 1 1 points). We 
use this spectrum estimate (P[ft)]) to compute the standardized 
and whitened residuals for each time series x, ■[/]: 

e ,[f]=iFFT(x,[co]/V^R) ■ (5) 

Here Xj [to] = FFT (xj [t]j is the Fourier transform of the original 
time series Xi[t) and iFFT(*) is the inverse Fourier transform 
of *. Finally, for each bootstrap replicate, we resample the 
values ei[t] with replacement and compute the surrogate data 

x i [t]=iFFT(e i [(0]- V^R) , (6) 

where e, [co] is the Fourier transform of the residuals e, [f] re- 
sampled with replacement. This last step ensures that the 
spectral characteristics (e.g., l// a behavior) of the original 
data are preserved in the surrogate data. 

We compute N s instances of these surrogate data, and for 
each instance we calculate the test statistic £y for each pair of 
nodes ; and j, (i.e., we calculate the maximum of the absolute 
value of the cross correlation between the surrogate data Jc; [t] 
and X/ [f ]). The N s instances of sn form a bootstrap distribution 
of maximum cross correlation values to which we compare sn 
observed in the original data and assign a p-value. 

Constructing the bootstrap distribution of Sy values for 
all node pairs is computationally expensive. If our network 
contains 100 nodes, then we would like to compute a boot- 
strap distribution (and test the significance) of each of the 
100 x 99/2 = 4950 values sy. If each bootstrap distribution 
requires N s = 10000 surrogates, a standard choice in the liter- 
ature, then we construct surrogate data and compute the cross 
correlation over 10 7 times. We reduce this expensive oper- 
ation in the following way: instead of computing bootstrap 
distributions for all electrode pairs, we compute the bootstrap 
distribution (with N s surrogates) for only a subset of node 
pairs. We then define the merged distribution as the combined 
distribution for the entire subset of node pairs. We use the 
merged distribution to test the significance of sy for all node 
pairs (even pairs not used in calculating the merged distribu- 
tion). Note, however, that in doing so we assume that the null 
distribution of sy is the same for all node pairs. 

C. Step 3: Control of the false detection rate 

To test the statistical significance of the values sy, we may 
apply either of the two methods described above (or even, 
practically speaking, the naive method as well). Networks 
of, say, 100 nodes will consist of 100 x 99/2 = 4950 sy val- 
ues, each with an associated /9-value. Clearly, multiple testing 
is an important concern. If we simply choose a standard p- 
value cutoff for assessing significance (such as p < 0.05), then 
we expect the number of network edges incorrectly declared 
present to scale proportionally (i.e., roughly 250 such edges, 
for an 0.05 cutoff). To control for this abundance of false pos- 
itives, we could define a stricter cutoff; for example, we could 
use the Bonferroni correction and divide the /?-value threshold 
by the number of node pairs (i.e., p < 0.05/4950 = 10~ 5 ). 



This conservative control of the familywise error rate (the 
probability of making one or more false discoveries) is likely 
too strict for data in which we expect relatively few significant 
edges a priori, i.e., for sparse networks. 

Instead, we employ the less conservative false discovery 
rate (FDR) to control for multiple testing. The FDR is de- 
fined as the expected proportion of erroneously rejected null 
hypotheses among the rejected ones ll49"[|50Tl . and various pro- 
cedures exist for controlling the FDR in practice ll32ll . Gener- 
ally speaking, the notion of FDR control guarantees that the 
expected proportion of falsely declared edges in our inferred 
networks is no more than a pre-specified fraction q € (0, 1). 
However, in order for this guarantee to hold, two assumptions 
must be true, namely that (i) statistical /^-values associated 
with each test are computed accurately, and (ii) tests are in- 
dependent. Of these assumptions, the first is critical, while 
the second is less so. That the second is less critical is im- 
portant in the context of network inference, since the various 
tests for declaring presence or absence of edges are clearly 
correlated, as they reuse the same time series. Additionally, if 
one wishes to address this dependency, there are extensions of 
the basic FDR procedure (e.g., see 113011 for useful discussion), 
although we do not pursue this here. On the other hand, in- 
accurate calculation of p-values is known to be disastrous to 
FDR principles. Our analyses presented below confirm this in 
the context of our network inference problem, and the major- 
ity of our efforts focus around this point, as we described in 
Section HUB] 

Here we implement the linear step-up FDR controlling 
procedure of Benjamini and Hochberg 114911 . which is com- 
puted as follows. First, order the m = N(N — l)/2 p-values 
Pi < P2 < ■•• < Pm- Then, choose a desired FDR level q. Fi- 
nally, compare each pi to the critical value q ■ i/m and find the 
maximum i (call it k) such that p^ < q ■ k/m (and therefore 
Pk+i > q • (k + 1 )/ m )- We reject the null hypothesis that time 
series x, -[t] and xj[t] are uncoupled for p\ < ■■■ < Pk- 

The choice of q determines a threshold p-value, for which 
we declare all features significant up through that thresh- 
old lIBTll . The value q represents an upper bound on the ex- 
pected proportion of false positives among all declared edges 
in our inferred network (i.e., among all node pairs for which 
was declared to be significant.) For example, if we fix 
q = 0.05 and find 100 significant values of \y, then we expect 
0.05 ■ 100 = 5 false positives (i.e., five false edges in the 100 
edge network). 



IV. RESULTS 

We analyze three data sets using the procedure defined in 
Section [HI] Two data sets we create with specific (known) 
structural topologies, to which we compare the functional 
topologies extracted through analysis of the dynamic data. 
The third we observe from a human epileptic subject undergo- 
ing invasive electrical monitoring of the cortex during seizure. 
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A. Pink noise data 

Many time series data produced by natural systems pos- 
sesses a l/f a power spectrum ll52l p3l Isill . To mimic this 
behavior, we create a nine-node network, first by generating 
500 ms (sampling interval 1 ms) of independent colored noise 
(a = 0.33) data w, [t] at each node ;. We then connect node i to 
j by adding pointwise to Wj [t] the data w, [t] scaled by a factor 
of 0.4. For example, we connect node #1 to #2 by adding to 
W2 [t] the time series 0.4 • w\ [t] for each time point t to create 
X2 [t] = W2 [t] + 0.4 • w\ [t], the time series associated with node 
#2. In Fig. 12 a) we illustrate the topology of the constructed 
network; a total of nine directed connections exist. 



(a) True (b) Naive (c) Extremum (d) Bootstrap 




Index Index Index 



FIG. 2: (Color online) For the colored noise network all three mea- 
sures perform equally well and detect the underlying structural topol- 
ogy. In the upper row, each circle indicates a node, (a) The true con- 
nectivity of the network. With colored (or shaded) lines we indicate 
directed connections between nodes; connections initiate from the 
red (dark) line end and terminate at the yellow (light) line end. (b-d) 
The p-values (lower) and corresponding functional network topolo- 
gies (upper) derived from the naive method (b), extremum method 
(c), and bootstrap method (d). The dark gray line in the lower fig- 
ures indicates the threshold for the linear step-up FDR procedure; 
we consider p-values below this line — in the unshaded region — 
significant. All three significance tests capture the functional topol- 
ogy equally well. 

Having established the network topology, we now attempt 
to recover it directly from the time series data. To do so, we 
apply our coupling measure (sy, the maximum of the absolute 
value of the cross correlation) pairwise to all m = 9 ■ 8/2 = 36 
electrode pairs in the network. We then test the significance 
of each sn and compute a p-value using the analytic and 
computational procedures defined above. We begin with the 
naive method, whose p-values we plot as asterisks in the 
lower portion of Fig. |2jb). Plotted in increasing order, these 
p-values range from ~ 10~ 5 to 0.3. We fix q = 0.10 and 
also plot in the lower portion of Fig. |2jb) the line of slope 
q/m = 0.10/36 = 0.0028 and zero intercept. Following the 
linear step-up FDR procedure, we reject the null hypothesis 
of no coupling for those (nine) electrode pairs with p-values 
below this line. We plot in the upper portion of Fig. [2fb) the 
(nine) "significant edges" corresponding to the significant p- 
values. Our confidence in this nine node network — derived 
from the time series data — is high; from the FDR procedure 



we expect 0.10 ■ 9 ~ 1 false positive edge (i.e., one spurious 
edge between uncoupled nodes). In this case, we find exact 
agreement between the known network topology (Fig. |2|a)) 
and the derived topology. We note that, for sake of clarity, 
we chose a simple coupling measure that does not determine 
edge direction. More sophisticated coupling measures that in- 
dicate edge direction may be employed following the general 
paradigm outlined above, as we discuss in SectionlVl 

In Figs. |2|c) and|2jd), we show the topology derived using 
the extremum and bootstrap methods, respectively. In both 
cases, we follow the linear step-up FDR procedure with q = 
0.10 to identify significant edges. For the extremum method 
(Fig. |2 C )) we identify eight significant edges, one less than 
expected. We compute our confidence in the network using 
the FDR procedure and anticipate 0.10 -8 ~ 1 false positive 
edge. 

To compute the frequency domain bootstrap, we first cal- 
culate the average power spectrum of all (nine) nodes. We 
then create a merged distribution using a subset of ten elec- 
trode pairs (of the possible 36) and N s = 10000 for each sur- 
rogate distribution. The resulting merged distribution con- 
tains 10 -N s = 10 5 points; therefore, the smallest p-value we 
can compute through this method is 10~ 5 . We find, in this 
case, six p-values at this detection limit. Using the bootstrap 
method (Fig. [2jd)), we identify ten network edges, one more 
than expected. We do expect 0.10 • 10 = 1 false positive edge 
in the network, although given only the time series data, we 
could not identify which of the ten edges is spurious. 

These simulation results suggest that all three measures of 
edge significance perform equally well. This is surprising, 
especially for the naive method in which we neither Fisher 
transform the maximal correlation values (to induce normal- 
ity), nor account for our choice of an extremum (the maxi- 
mum of the absolute value of the cross correlation). The naive 
method succeeds, in this case, because the two omissions ap- 
pear to balance. Omitting the Fisher transformation increases 
the p-values we observe, while utilizing the normal distribu- 
tion with zero mean — not the extremum distribution — de- 
creases the p-values. One omission compensates the other so 
that, in this case, the resulting p-values are approximately cor- 
rect. Unfortunately, we cannot rely on this delicate balance to 
always succeed as we illustrate in the next example. 



B. Simulated neural data 

In the previous model, we simulated colored noise activity 
possessing a 1 / f a falloff of the power spectrum. We now con- 
sider a more realistic model of interacting neural populations. 
We provide a brief description of the model here; more details 
may be found in Appendix 2. The model consists of 1000 neu- 
rons divided into twenty groups of 50 cells. Within each group 
we include strong connections (excitatory synapses) between 
randomly chosen neurons; activity initiated by a few neurons 
in a group quickly spreads to the other neurons of the same 
group. Between cell groups, we establish only weak (excita- 
tory synaptic) connections joining individual neurons of spe- 
cific cell groups. We illustrate the topology of these weak 
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connections between cell groups in Fig. EJa). In this figure, 
each gray circle represents a cell group (of 50 strongly con- 
nected neurons) and lines represent connections between cell 
groups. With this connectivity in place, we simulate the neu- 
ral dynamics and compute the average population activity of 
each group. We then employ the general paradigm described 
above to the resulting neural activities and compare the mea- 
sured functional connectivity (i.e., the pattern of connectivity 
inferred from the simulated neural dynamics) to the known 
structural connectivity between nodes shown in Fig.[3ja). The 
results, as we show below, depend upon the significance test 
we employ. 



(a) True (b) Naive (c) Extremum (d) Bootstrap 




FIG. 3: (Color online) For the simulated neural data the choice of sta- 
tistical test is vital to construct an appropriate network, (a) The data 
consist of twenty cell groups (gray circles) and 22 connections be- 
tween cell groups. Directed connections proceed from the red (dark) 
to yellow (light) end of each line, (b-d) The functional networks de- 
duced. The naive method (b) identifies no significant p-values; with 
q = 0. 10 in the linear step-up FDR procedure, none of the p- values 
lie below the (dark gray) line (q/m) ■ i. The extremum method (c) 
identifies 17 significant edges (of which we expect two are false pos- 
itives); 14 match the structural network in (a). The bootstrap method 
(d) detects 18 edges, of which we expect 2 false positives. This pro- 
cedure detects 15 (of the 22) true edges. 

We apply the coupling measure pairwise to all m = 20 • 
19/2 = 190 possible group pairs in the network and test the 
significance of each result by computing a /9-value using one 
of the three procedures defined above. We begin with the 
naive method, whose /^-values we plot as asterisks in the lower 
portion of Fig.|3jb). With q = 0.10 in the linear step-up FDR 
procedure, we find no significant values of maximal cross cor- 
relation; none of the /^-values lie below the line (q/m) ■ i. 
The resulting (trivial) network — shown in the upper por- 
tion of Fig. [3b) — contains no edges. The other two signifi- 
cance tests produce nontrivial networks. Using the extremum 
method and linear step-up FDR procedure (with q = 0.10) we 
identify 17 significant edges. The resulting network, shown 
in Fig. |3]c), correctly identifies 14 edges and possesses three 
erroneous edges (i.e., edges we identify in the functional net- 
work that do not exist in the structural network). We expect 
from the FDR procedure q ■ 17 ~ 2 false positives, in approxi- 
mate agreement with the three erroneous edges observed. Fi- 
nally, we show in Fig. |3jd) the /^-values and network deter- 
mined using the bootstrap method. In this case, we detect 18 



edges (and expect 2 false positives). This procedure detects 15 
(of the 22) true structural edges and produces three erroneous 
edges, again in approximate agreement with the number of 
false positives expected. 

In all three cases, the functional topology derived from 
the mean dynamics fails to capture exactly the true structural 
topology of the network. The naive method detects no signif- 
icant edges and performs most poorly. This is not surprising; 
we expect that the un-normalized /^-values and incorrect dis- 
tribution of maximal correlation values will compromise the 
naive method. The extremum and bootstrap methods produce 
similar functional networks that approximate the true struc- 
tural network. Although both measures make mistakes, the 
FDR procedure provides an estimate for the number of er- 
roneous edges to expect. We conclude that, for these simu- 
lated data, the extremum and bootstrap methods outperform 
the naive method and qualitatively reproduce many (but not 
all) of the network edges. 



C. Human ECoG data 

In the previous two examples, we applied the coupling anal- 
ysis to networks with known structural topology. This al- 
lowed us to compare the derived functional topology with the 
true structural topology and determine each method's perfor- 
mance. As a last illustration of the methods, we consider volt- 
age activity recorded directly from the cortical surface (elec- 
trocorticogram or ECoG data) of an epileptic human subject 
for clinical purposes (Appendix 3). We focus on a short inter- 
val (1 s) of data recorded from 97 electrodes while the subject 
experienced a seizure. We apply all three methods to the data 
and compare the resulting (functional) networks. In this case, 
the structural connectivity is unknown. We find that, as be- 
fore, the extremum and bootstrap methods produce consistent 
results. 

We show the deduced functional networks in Figs. HJb- 
d). In each case, we test the significance of m = 97 • 96/2 = 
4656 maximal cross correlation values, and use a linear step- 
up FDR procedure with q = 0.05 to define significant p- 
values. For the naive method (Fig. SJb)) we find no signifi- 
cant /^-values and the corresponding trivial network contains 
no edges. We note that the node locations in Fig. [4] do not 
correspond to their physical locations on the human cortex. 
Instead, we simply arrange the nodes in a circle. 

From the extremum and bootstrap methods we create simi- 
lar networks. For the former, we identify 162 significant edges 
(of which we expect 9 false positives) as drawn in Fig. |4jc). 
For the latter, we select 500 electrode pairs (of the possible 
4656 pairs) to compute surrogate distributions, each distribu- 
tion containing N s = 10000 realizations. The smallest p-value 
detectable in the resulting merged distribution is 2 x 10~ 7 . Us- 
ing this method we find the 187 significant edges drawn in Fig. 
HJd), of which we expect 10 false positives. 

Comparing the functional networks deduced from the ex- 
tremum and bootstrap methods, we find that the two are sim- 
ilar. Moreover, we show in Fig. |4fa) a fourth functional net- 
work constructed using a simple threshold procedure; we in- 
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elude edges only between those node pairs with su > 0.75. 
Remarkably, all three networks are qualitatively similar al- 
though we use different techniques to construct each network. 
Of course the simple threshold network does not indicate our 
confidence in the network: how many edges in Fig. |4{a) are 
false positives? In addition, we note that the bootstrap method 
is computationally expensive; constructing the surrogate dis- 
tribution requires approximately 90 minutes on a 2 GHz Core 
Duo processor and therefore at least 45 hours to construct the 
networks for 30 s of seizing activity. The extremum method, 
designed for our particular choice of coupling measure, iden- 
tifies a network similar to the bootstrap method in a computa- 
tionally efficient way. 



(a) Simple (b) Naive (c) Extremum (d) Bootstrap 




FIG. 4: Functional networks constructed from 1 s of ECoG data 
recorded at 97 electrodes during a seizure depend upon the statis- 
tical test we perform, (a) A simple threshold network with edges 
(black lines) drawn between nodes pairs exhibiting sufficient func- 
tional coupling, Sij > 0.75. (b-d, lower) The 4656 p-values calcu- 
lated from the naive method (b), extremum method (c), and boot- 
strap method (d). For each method, we fix q = 0.05 in the FDR 
procedure, (b-d, upper) The corresponding functional networks. The 
naive method (b) detects no significant edges and the corresponding 
network is trivial. The network created from the extremum method 
(c) contains 162 edges, and from the bootstrap method (d) 187 edges. 



We follow the procedure described above to analyze these 
"shuffled" data. We compute the maximal cross correlation 
for each electrode pair, and show the corresponding p-values 
and functional networks in Fig.|5jb-c). With q = 0.05, we find 
no significant coupling using the naive or extremum methods. 
We do detect 2 significant edges with the bootstrap method 
(of which we expect 1 false positive). These significant edges 
match those determined using a simple threshold procedure 
(sij > 0.75) whose network we show in Fig. EJa). We con- 
clude that the three significance tests behave as expected for 
the shuffled data; if we disrupt the coupling in the data, we 
expect trivial functional networks. 



(a) Simple (b) Naive (c) Extremum (d) Bootstrap 




Index index index 



FIG. 5: By shuffling the ECoG data, we eliminate coupling between 
the time series and detect no (or few) edges, (a) A simple threshold 
network with edges drawn between node pairs with sufficient func- 
tional coupling (sij > 0.75) detects two edges located at the left of 
the network, (b-c) The p-values (lower) and corresponding networks 
(upper) derived from the (b) naive, (c) extremum, and (d) bootstrap 
methods. Only the latter detects two edges (of which we expect 1 
false positive). 



V. DISCUSSION 



D. Human ECoG data: shuffled 

For the human ECoG data, we do not know the structural 
network (i.e., we do not know the topology of chemical and 
electrical connections between neurons in these cortical re- 
gions). Therefore, we cannot validate the functional networks 
shown in Fig. |4]by comparison with anatomical connections. 
However, we can manipulate the ECoG data to disrupt func- 
tional connections and verify that our significance tests detect 
no coupling. To do so, we create a new data set: we assign 
to each electrode 1 s of data chosen at random from a 120 s 
interval that includes 60 s of pre-seizure and 60 s of seizure ac- 
tivity. For example, electrode #1 may contain ECoG data from 
t = [8 .2, 9.2] , electrode #2 data from t = [97 .0, 98 .0] , electrode 
#3 from f = [1 10.4, 1 1 1 .4], and so on. With the data chosen in 
this way, we expect only weak associations between electrode 
pairs. 



Our increased ability to collect multivariate spatiotempo- 
ral data (e.g., from high density electrode arrays) necessitates 
the construction and analysis of complex functional networks. 
In this manuscript, we adopted a statistical hypothesis test- 
ing paradigm for constructing such functional networks. This 
paradigm involved three steps: 1) choice of an association 
measure, 2) definition of a significance test, and 3) accounting 
for multiple significance tests. Although the paradigm itself is 
quite general, the details accompanying each step are problem 
specific. 

Here we developed this general paradigm for multivariate 
time series data. For the association measure we chose the 
maximum of the absolute value of the cross correlation. We 
defined two approaches to significance testing (one analytic 
and the other computational), and employed a linear step-up 
FDR procedure to account for multiple tests. Applying these 
techniques to three data sets, we showed that the choice of sig- 
nificance test was critical. Without accurate /^-values for each 
network edge, we lack confidence in the resulting network. 
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The general paradigm outlined in Section [XT] applies to any 
choice of association measure. In this work we focused on 
this simple cross correlation measure for two reasons. First, 
the measure is computationally efficient. Second, analytic ex- 
pressions exist (or can be derived) to test the significance of 
each cross correlation result. More appropriate coupling mea- 
sures exist ll35tl that may perhaps improve the network results 
we present here. In particular, measures that distinguish direct 
from indirect interactions and incorporate the flow of informa- 
tion ll55l l56ll would be of use. However, choosing a more so- 
phisticated association measure does not guarantee more ac- 
curate functional networks. The coupling measure must also 
include an accurate significance test; without precise p-values 
for each network edge, we weaken our measures of network 
confidence. 

Researchers in various other contexts have followed a sim- 
ilar strategy of associating p-values with each network edg e 
and accounting for multiple significance tests (e.g., I57ll58l0 . 
Our numerical results illustrate how the choice of an appropri- 
ate significance test associated with a specific coupling mea- 
sure is critical. That a measure possesses a significance test 
does not guarantee accurate p-values; typically significance 
tests make specific assumptions about the data. For exam- 
ple, we found that the naive method — although perhaps in- 
tuitively appealing — was inappropriate because we did not 
account for taking the maximum of the absolute value of the 
cross correlation, and thus produced inaccurate p-values and 
inaccurate networks. Therefore, we utilized two additional, 
complimentary measures. By testing the paradigm on simu- 
lated data with known physical connectivity we deduced ap- 
propriate significance tests for the association measure imple- 
mented here. 

We note that trivial networks (e.g., networks without edges 
as in Fig. [5]l rarely appear in practice. Upon finding a trivial 
network, a common response is to adjust the network thresh- 
old to include more edges, perhaps until the network becomes 
strongly connected. To follow a similar strategy here we in- 
crease the value of q in the linear step-up FDR procedure. If, 
for example, we set q — 0.5 (instead of q = 0.05) we may de- 
tect new significant network edges. But by increasing q we 
decrease our confidence in the network; with q = 0.5, we ex- 
pect half of the network edges declared significant are false 
positives. Thus, through our choice of q, we balance the num- 
ber of edges detected with our confidence in the network. 

The typical approach to construct functional networks from 
multivariate time series data involves thresholding an associ- 
ation measure. For example, we may define edges between 
nodes whose maximal cross correlation exceeds 0.75, as in 
ll22ll . This procedure for constructing a network suffers from 
numerous inadequacies. First, we lack a measure of confi- 
dence in the resulting network. With this choice of 0.75 as 
threshold how many spurious edges do we expect, and does 
this number change as we vary the threshold? Second, we 
expect the choice of threshold may depend on the particu- 
lar instance of data observed. For example, in constructing 
functional networks of ECoG data recorded during seizure, 
the threshold may vary from patient to patient, depending on 
mechanisms intrinsic to each individual. Finally, a more ro- 



bust approach to constructing functional networks must prop- 
agate error in the association measure to uncertainty in net- 
work measures (e.g., to uncertainty in measures of degree or 
betweenness). 

We propose that choosing a threshold value of q, rather 
than a threshold value of an association measure, constitutes a 
more rigorous procedure for establishing functional networks. 
By choosing the threshold through the use of formal statistical 
hypothesis tests, we create functional networks with specified 
levels of network uncertainty that may be calibrated across 
a population of multivariate data. In the future, we will use 
this approach to study how uncertainty in the association mea- 
sure affects uncertainty in network characteristics, and how 
to adopt these procedures for weighted (rather than binary) 
networks. Combined with biophysical models, robust tech- 
niques to create functional networks will perhaps illuminate 
the mechanisms that produce the observed activity and, when 
necessary, suggest how to alter this activity. 
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VL APPENDIX 1: DERIVATION OF ©. 

Suppose that Z\ , . . . ,Z„ are independent and identically dis- 
tributed normal random variables, with mean and variance 
1. Define M„ = max, (Z, ) and m„ = min, (Z, ). Then 

Pr (max |Z,| < zj = Pr (M„ < z , m„ > -z) 

w Pr(M„ <z)Pr(m„ > -z) 
= [Pr(M„<z)] 2 . 

The approximate equality in the second line follows by 
asymptotic independence of the max and min (e.g., Theorem 
1.8.2 of I59I0 . while the equality in the last line follows by 
symmetry of the normal distribution. Now by, for example, 
Theorem 1.5.3 of l59ll we have 

Pr (a n (M„-b n ) <z)«exp(-e- z ) , (7) 

with equality holding asymptotically in n. As a result, we 
obtain the expression in (0]i in the ideal case that the standard- 
ized statistic zfj derives from cross correlations CC,y [x] that are 
independent. Although these cross correlations will clearly 
be dependent, the approximation (|4]i nevertheless can be ex- 
pected to hold fairly generally, as the basic limiting extreme 
value distribution used here is quite robust for sequences of 
normal random variables under a range of dependency condi- 
tions, for both stationary and even non-stationary cases. See, 
for example, Chapters 4 and 6 of JUt]. 
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VII. APPENDIX 2: NEURAL MODEL 

We model the dynamics of each neuron with two ordinary 
differential equations, one to represent the membrane volt- 
age, and the other a membrane recovery variable l6(ill . We 
choose the model parameters so that each neuron generates 
regular spiking activity (i.e., a = 0.02, b = 0.2, c = —65.0, 
d = 8.0 in [60]). We then connect the neurons with excita- 
tory synaptic connections to establish two connectivity pat- 
terns: strong-local connectivity and weak-distant connectiv- 
ity. In both cases, we divide the neurons into groups of 50 cells 
numbered sequentially (i.e., group #1 contains cells {1 — 50}, 
group #2 cells {51 - 100}, group #3 cells {101 - 150}, and 
so on.) Within each local group of cells, we create 1200 di- 
rected excitatory synapses (of the possible 50 x 49 = 2450 di- 
rected pairs with no self synapses). Each synapse is assigned 
a uniform random conduction delay between and 10 ms and 
synaptic strength chosen uniformly between and 15. These 
synapses establish the strong-local connectivity within a cell 
group and define the twenty cell groups in the network; see 
Fig.! 

We also create weaker synaptic connections between the lo- 
cal cell populations. To do so, we select two groups (e.g., #1 
and #8) and create 550 excitatory synapses from neurons in 
one group (e.g., #1) to neurons in another (e.g., #8). These 
"distant" synapses are weaker than the local connections; we 
assign the synaptic strengths smaller random values (chosen 
uniformly between and 5) and uniform random conduction 
delay between and 10 ms. We illustrate the 22 distant con- 
nections between the twenty cell groups in Fig. 0a). Each 
gray circle represents a local cell group (i.e., a subset of 50 
neurons). The colored (shaded) lines represent the distant 
synaptic connections between groups. In addition to the local 
and distant synaptic inputs, we also include strong synaptic in- 
put (from the thalamus, say) to one randomly chosen neuron 
each millisecond, causing this neuron to generate an action 
potential. We follow the algorithm in ll60ll to simulate the neu- 
ral population for 5000 steps (or 5 s) with a sampling interval 
of 1 ms. The model is similar to recent simulations ll6lll62ll . 
except that we introduce here connectivity with a particular 
structural topology. 

In the human ECoG recordings described in Section |TV C| 
we observe the dynamics of postsynaptic potentials produced 
by large neural populations, not the spiking activity of indi- 
vidual neurons 11631 l64ll . To mimic these population dynam- 
ics, we construct the mean activity of the local neural groups 
in the following way. First, we define I[t] as the total current 
input to each neuron at time t . In the model equations we sim- 
ulate here, these current inputs change instantaneously l60ll . 
In reality, current inputs follow the opening and closing of 
channels and evolve more slowly l65ll . To approximate these 
slow dynamics at postsynaptic neuron j, we use the following 
equation: 

Sj=Ij[t]{l-Sj)-sj/20, (8) 

where Sj represents the state of the synapse at neuron j, and 
Ij [t] represents the total excitatory input current to the neuron 
j at time t. We note that, for simplicity, we approximate the 



total synaptic input to neuron j as a single synapse with dy- 
namics driven by Ij[t], the activity of all neurons presynaptic 
to j. When Ij [t] is large, excitatory current enters neuron j and 
sj — > 1. When Ij[t] is small, sj — ► with a decay time constant 
of 20 ms, and neuron j approaches its resting potential. We 
define the mean activity of a neural group as the average of Sj 
over the local (fifty) cell population. For example, we com- 
pute the mean activity of population #1 as the average value 
of sj for neurons j = {1,2, ...,50}. We only use Sj to define 
the mean population activity; these synaptic dynamics do not 
impact the voltage dynamics of the model neurons. 

With the population activity defined in this way (i.e., as the 
mean synaptic dynamics within a neural group), we simulate 
5 s of neural dynamics and record the average activity of each 
group. We then scale the group activity to have zero mean 
and unit variance and add Gaussian noise (zero mean and 
0.55 variance) to each sample of each time series. Finally, 
we downsample the group traces by a factor of five, reducing 
the sampling rate (from 1000 Hz to 200 Hz) to decease subse- 
quent computational time; we show examples of the resulting 
time series data in Fig. [6] 

rfi=TTli^r>T^ ... ^krn 
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FIG. 6: (Color online) A cartoon representation of the network used 
to simulate the neural data. Each group of neurons contains 50 cells, 
represented as filled triangles. Within a group, we include many, 
strong synaptic connections (terminating at large circles, orange). 
Between groups, we include few, weak synaptic connections (ter- 
minating at small circles, blue). We show examples of the average 
activity of each group (from which we construct the functional net- 
works) in the lower portion of the figure. 



VIII. APPENDIX 3: HUMAN SUBJECT DATA 

The ECoG data were recorded from a 37 year old male pa- 
tient with medically refractory epilepsy whose seizures began 
at age 3. Following the failure of seven antiseizure medica- 
tions and a vagal nerve stimulator, and upon recommendation 
of his clinical team consisting of epileptologists and neurosur- 
geons, the decision was made to pursue resection of the tissue 
from which the seizures arose. To this end, subdural grids of 
electrodes were implanted. The goal of this procedure was 
to identify the epileptogenic zone — the region of the brain 
producing recurrent seizures — and surgically remove it j66ll . 

The ECoG recordings consisted of 100 electrodes placed 
directly on the cortical surface (92 electrodes over the left 
frontal and temporal lobes) and within deep brain regions (8 
electrodes within the temporal lobe). Following electrode im- 
plantation, the subject was admitted to a specialized moni- 
toring unit and data recorded continuously at 500 Hz for ten 
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days. During this time, four seizures were observed; to illus- 
trate the measures, we analyze only the second seizure here. 
Analysis of the collected data was approved through the Part- 
ners Health Care Human Research Committee and the Charles 
River Campus Institutional Review Board. 




100 ms 



FIG. 7: (Color online) Examples of the ECoG data recorded from 
the human subject during seizure. In the upper trace we show 160 
s of data recorded from a single electrode; the seizure begins at the 
transition from low amplitude to high amplitude fluctuations denoted 
by the arrow. We indicate the one second interval analyzed with the 
vertical red line. The lower ten traces illustrate the voltage activ- 
ity recorded from ten (of the 97) electrodes during this one second 
interval; we apply the coupling measure to all pairs of these data. 



We apply our coupling analysis to simultaneous record- 
ings from 97 electrodes; three electrodes — suffering from 
extreme artifacts — were discarded. Before beginning the 
coupling analysis, we process the ECoG data in the following 
way. First, we lowpass filter the data (two-way least-squares 
FIR filtering) below 55 Hz to isolate the low frequency com- 
ponents. We therefore ignore higher frequency activity that 
may delineate seizure onset 11671 16811 and instead focus on the 
high amplitude, low frequencyoscillations that characterize 
unequivocal clinical seizures Il69ll . We then choose a 1 sec- 
ond interval of the ECoG data during the seizure. We chose 
this short interval to balance two competing needs: station- 
arity and sufficient data. For the former, we must choose an 
interval in which the voltage dynamics at each electrode re- 
main approximately consistent (i.e., exhibit oscillations of the 
same approximate character). For the latter, we must choose 
an interval that contains enough data to calculate the coupling 
measure (e.g., a 50 ms interval would fail to capture some 
slow oscillations characteristic of a seizure). Finally, we nor- 
malize the data from each electrode within the 1 second inter- 
val to have zero mean and unit variance. We show examples 
of the ECoG data employed in the analysis in Fig. [7] 
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